Analysing continuous proportions in ecology and evolution: A practical introduction to beta and Dirichlet regression

3.1 Analysing proportions with two continuous categories

3.1.2 Beta regression

Beta regression is a technique that has been proposed for modelling of data for which the observations are limited to the open interval (0, 1) (Ferrari & Cribari‐Neto, 2004; Smithson & Verkuilen, 2006). Some recent examples of beta regression in ecological contexts include analyses of: the contribution of food derived from different energy pathways (Child & Moore, 2015; Fukumori, Yoshizaki, Takamura, & Kadoya, 2016); forest simulation output in terms of proportions of tree basal area belonging to focal tree species (Ameztegui, Coll, & Messier, 2015); and the degree of leaf damage due to a leaf pathogen (Busby et al., 2013). Although we present beta regression in this review as a method for analysing continuous‐based proportions, there are also many examples for using this techniques to analyse data from derived indices that are bound between 0 and 1 (e.g. an evenness index in Nogueira, González‐Troncoso, and Tolimieri (2016); or a straightness index in Shimada et al. (2016)). In addition, it can be applied to variables that are constrained to the interval a and b as they can be rescaled to [0, 1] through (y − a)/(b − a).

Beta regression consists of the same three components as generalized linear models (GLMs) (Bolker et al., 2009; McCullagh & Nelder, 1989), and those familiar with GLM will recognize the most important aspects of beta regression (the distinction between the two arises from the non‐orthogonality of the model parameters, see below). Here, we briefly review these three elements: the random component (the beta distribution and its implied mean–variance relationship), the systematic component (the linear predictor) and the link function (specifying the link between the random and systematic component). We refer readers to Ferrari and Cribari‐Neto (2004) for a more comprehensive explanation.

Beta regression begins with the assumption that the data‐generating process can reasonably be modelled by a beta probability distribution (Balakrishnan & Nevzorov, 2003). The beta distribution is a member of the exponential family (Kieschnick & McCullough, 2003), and is defined by two parameters for values on the open interval (0, 1). Two parameterizations for the beta distribution are available, but the mean‐precision parameterization, with µ (for the expected value) and ϕ (as a measure of ‘precision’, or the inverse of dispersion), is most commonly used in the context of beta regression (see Box 1). The variance can be related to the mean (µ) by and is therefore proportional to the variance of the binomial distribution for one trial, µ(1 − µ), by a factor of .

Depending on the choice of values for the two parameters a large range of shapes can be obtained including symmetrical, skewed, uniform, roughly bell‐shaped and bimodal. This flexibility, combined with the limitation to values between 0 and 1, make the beta distribution a particularly useful model for continuous proportional data. In addition, fitting a beta distribution gives increasingly less biased estimates of the mean compared to transformation‐based approaches when observations get closer to zero and one and/or their variance is large (see Appendix S1).

BOX 1. Mathematical details of beta and Dirichlet regression

Beta regression

Definitions of the beta distribution usually employ a parameterization using the symbols α and β such that the probability density function for a beta‐distributed response variable y is given by the following:

where

And Γ(.) is the gamma function.

The corresponding expectation and variance of the distribution are given by the following:

When using the beta distribution for modelling data, it is usually more convenient to use an alternative parameterization with μ and ϕ, such that the expectation of the distribution is simply E[y] = μ, and the variance is given by the following:

In beta regression, the conditional model for the mean μ of the response given covariates X is usually assumed to be linear on the logit transformed scale:

where η is known as the linear predictor, β is a vector of parameters to be estimated and X is the design matrix of covariate values.

To obtain estimates from a beta regression fit that are interpretable on the scale of observations (0, 1) the values from the linear predictor therefore need to be back‐transformed with the inverse logit function:

As discussed in the text, ϕ can either be estimated as a single value for all observations, or modelled as a function of covariates with design matrix Z, corresponding regression parameters γ, and linear predictor ζ – in which case a log link is appropriate:

Dirichlet regression

The traditional parameterization results in the following probability density function for a vector valued p.

where x and α are both vectors (of observations and model parameters respectively) of length C, and B() is the multinomial Beta function.

In this parameterization, there is a parameter α_c for each of the C components. If we define α₀ as the sum of all elements of α then the expected value of any given component x_c is given by , with the associated variance:

The value of α₀ is therefore interpretable as a precision parameter.

An alternative parameterization is realized by representing the expected values of each of the components as a vector μ where each μ_c is between 0 and 1, and the sum of all elements of μ is 1. Additionally, we define a precision parameter ϕ. Conversion between the two paramaterizations is therefore possible with α_c = μ_cϕ and α₀ = ϕ, leading to the follow expressions for expected value and variance of each component:

The main contrast with the traditional parameterization is therefore that the group means and the precision parameter are explicitly modelled rather than indirectly in terms of the values of α.

In the case of the traditional paramaterization, a regression model should be fitted for each of the C values of α_c. An appropriate link function for the corresponding regression model is the log‐function.

Under the alternative paramaterization, both the expected values μ and the precision parameter ϕ are modelled as a function of covariates. One of the c components is defined as the base category b (all regression coefficients = 0). The sum constraint implies that the appropriate link function for the regression models for μ is the multinomial logit function. In this way, given a separate linear predictor for each component η_c = X_cβ_c, these can be converted into the corresponding values of μ with the expressions:

Note that the expected value of the baseline component is implictly modelled through the models for the other C − 1 components. As for beta regression, a model for ϕ is typically specified with the log–link function.

Once a beta distribution has been chosen, the next step is the specification of the systematic component of the model relating the expected value of the response variable to one or more (continuous or categorical) predictor variables. In the familiar case of ordinary linear regression this dependence is specified through the regression equation µ = E[Y|X] = β₀ + β₁X₁ + ⋯ + β_pX_p where β are regression parameters to be estimated based on observed values of y and corresponding matrix X of p predictors. Given that such a linear predictor function can potentially range between −∞ and +∞ it cannot be used to specify the mean for distributions, such as the beta, that are restricted to a particular interval. For this reason, a link function must be specified to convert between the linear predictor model and the conditional mean on the scale of observations (Zuur et al., 2009). The model relating the values of the covariates, and the expected value of the response therefore becomes g(µ) = β₀ + β₁X₁ +⋯ + β_pX_p, implying that µ = g⁻¹(β₀ + β₁X₁ + ⋯ + β_pX_p), where g(.) and g⁻¹(.) are an appropriate link function and its inverse, respectively. In the case of beta regression, several link functions are potentially applicable, with the logit function (see Box 1) being the most common choice. The inverse of this function ensures that any value from the linear predictor will fall between 0 and 1. Appendix S2 presents a simple example of beta regression to make these concepts more concrete.

Estimates of β and ϕ that lead to a model that best fit the observed data are obtained by maximum likelihood estimation. The probability model (including covariates) and data are combined to define a likelihood function (see Bolker, 2008, for an accessible treatment), for which the maximum is obtain by numerical optimization, leading to maximum likelihood estimates for the parameters of the mean model (the βs) as well as for the precision parameter ϕ. The estimation procedure also produces values for the standard error of each parameter, allowing the application of the usual inference tools such as significance testing and confidence intervals. Here, it should be noted that for beta regression, the parameters β and ϕ are not orthogonal (Ferrari & Cribari‐Neto, 2004), complicating estimation for β when ϕ is not known, or incorrectly specified. This non‐orthogonality between the mean model parameters and the error parameter is also the reason why beta regression with unknown ϕ is not strictly a GLM (Huang & Rathouz, 2017). Moreover, two important properties of the beta distribution, not shared by common GLMs, are that the maximum likelihood estimate for µ can be different from the corresponding sample mean, in particular for small sample sizes; and that changes in the way that ϕ is modelled (see Section 4.1 below) can have consequences for the maximum likelihood estimates of µ. This can lead to bias in estimation and inference (see Section 5.1 and Case Study 1, below).

Although best suited for continuous proportions, beta regression has also been successfully employed to analyse count‐based proportions in cases where the numbers of observed units is large (e.g Bennett, Nimmo, & Radford, 2014; Briand, Schwilk, Gauthier, & Bergeron, 2015). However, a few cautionary remarks should be made. The standard error of a sample proportion decreases with the number of trials, n, according to . Thus, with a small number of trials random selection error of the trials may be important. In addition, when employing beta regression to count‐based proportions one loses information on the number of trials that was used to calculate the proportion. Effectively, each proportion is given equal weight, which can be problematic if the number of trials varies across samples. Two potential alternatives in this case would be to apply ‘beta‐binomial regression’ models (Skellam, 1948), or the use of an ‘observation‐level random intercept’ (Harrison, 2015). These methods are particularly suited to situations where count proportions show more variance than can be modelled by binomial regression, i.e. when count proportions are overdispersed. Overdispersion is more likely to occur when observations are grouped in space or time or when W becomes very large. In a beta‐binomial regression model, the probability of success for a given level of the covariate is not fixed, but comes from a beta distribution. In the binomial model with observation‐level random intercept each observation gets a random intercept. These random intercepts come from a Gaussian distribution at logit scale. A simulation study showed that in general a beta‐binomial is preferred and lead to least biased estimates (Harrison, 2015).

3.2 Analysing compositions with >2 continuous categories

4.1 Independent model for precision ϕ in beta and Dirichlet regression

In many situations, it is useful to allow for the possibility that precision varies across observations as a function of one or more covariates. Such a variable ϕ model may be appropriate when, for example, certain treatment levels display more variance at a given predicted mean (Case Study 1), or if there is systematic increase or decrease in variance as one of the covariates changes. It also allows for modelling situations where the variance of a response changes without any systematic change in the mean value. An additional use is for cases when differences in variability are a result of the sampling or experimental design, i.e. in cases where the observation interval cannot be fixed a priori. For example, in animal behavioural studies, the observation interval may be determined by how long the subject is within sight. In such cases, the precision of the observation can be expected to partly depend on the scale of the observation unit. Intuitively, an observation of 50% of time spent on a certain activity has a different evidential value when it is based on a 20 s total observation time, when compared to the same observed proportion derived from a 5 min observation interval. In such a situation, it may be appropriate to incorporate information about the length (or size) of the observation unit by including it in a model for the precision term ϕ.

Fitting beta and Dirichlet regression models with variable ϕ parameters is accomplished by extending the model and associated likelihood function to include dependence of ϕ on pre‐specified covariates (Simas, Barreto‐Souza, & Rocha, 2010; Smithson & Verkuilen, 2006). Since ϕ is defined to be always positive, a link function (usually the log function) is used to relate the continuous linear predictor to the ϕ scale. The fitting procedure will return maximum likelihood estimates and associated standard errors for the parameters of the model for ϕ which can be interpreted in the same way as the parameters for µ, after taking into account the different link functions.

4.2 Hierarchical data structures

Many experimental designs in ecology and evolution lead to observations that are grouped in some way. Common examples are multiple observations within an experimental plot or repeated observations of the same experimental unit through time. Such experimental designs lead to non‐independence of observations, violating assumptions of most statistical techniques and potentially leading to incorrect inference. Mixed effect models can account for non‐independence of observations and can be implemented as extensions to Generalized Linear Models (Bolker, 2015; Zuur et al., 2009) and have been developed for beta distributed variables as well (Brooks et al., 2017). Case study 1 includes a section on the fitting and interpretation of such a mixed‐effects model in the context of beta regression. At the time of writing, mixed effect models for Dirichlet regression have not yet been implemented in standardized software to the best of our knowledge, although see Regular et al. (2014) for an example of a customized analysis using Markov Chain Monte Carlo techniques to analyse seabird time budgets.

5.1 Bias in estimation

Beta and Dirichlet regression use maximum likelihood methods for parameter estimation. However, in many cases, the methods are known to lead to bias – a systematic deviation of an estimated parameter value from the true value – particularly in cases with small sample size (Firth, 1993). Biased estimators can lead to erroneous inference, and methods for bias–reduction and bias–correction are an active research area in statistics and in beta regression in particular (Grün et al., 2012; Kosmidis, 2014; Kosmidis & Firth, 2010; Ospina, Cribari‐Neto, & Vasconcellos, 2006; Simas et al., 2010). In particular, the precision parameter in beta regression models is prone to overestimation bias (Kosmidis & Firth, 2010) which leads directly to underestimation of the width of confidence intervals for other model parameters. Two main types of solutions are available to reduce bias: bias correction and bias reduction (Firth, 1993). Bias correction methods correct for bias in a separate step following maximum likelihood estimation, while bias reduction methods modify the maximum likelihood estimation procedure such that the resulting estimator is less biased. See Appendix S3 accompanying Case study 1 for a demonstration of bias correction and bias reduction and a bootstrap‐based technique for assessing the degree of bias for any given data‐model combination (Kosmidis, 2014).

5.2 Dealing with 0 and 1 in observations

The proportions modelled by the beta and Dirichlet distributions are defined on the interval (0, 1). However, for some combinations of parameters, the probability density function is zero or infinity at the boundaries, which precludes a meaningful calculation of the likelihood. Zeros and ones may occur in the data because the true (non‐zero) value is below the detection limit of the measurement device or method. Conversely, true zeros may occur when a given category is absent from the sample e.g. no vegetation is present in a sampled quadrat. Regardless of their source, observations of zero or one will lead to a failure of the fitting algorithms in both beta and Dirichlet regression. Note that this problem is not unique to these two techniques – both logit or additive log‐ratio transformations cannot be applied to datasets containing zeros and ones.

Here, we focus on the case of observations of zero, but the same advice and techniques can be applied to datasets containing observations of 1, or both. Several solutions are available depending on the origin of the zero. When the zero arises because of the detection limit of the observation method, a simple workaround is to replace all observed zeros with a small term prior to computing p for each sample, taking care to include ε in the new demoninator. ε can be chosen equal to the detection limit or to the smallest non‐zero observation. Warton and Hui (2011) recommend exploring the sensitivity of results to the value of ε.

Alternatively, the data can be transformed according to the following equation:

(1)

with p being the proportion of a category, n the total number of observations in the dataset, and C the number of categories (Maier, 2014; Smithson & Verkuilen, 2006).

Note that for a fixed value of n this is a linear transformation, and does not lead to issues arising from Jensen's inequality.

In case of exact zeros, as an alternative to the transformations above, it is possible to apply the so‐called zero‐inflated beta regression (Fang & Kong, 2015; Ospina & Ferrari, 2012), perhaps more properly referred to as zero‐augmented beta regression (Wright, Irvine, Warren, & Barnett, 2017), given the independence of the two processes. This type of regression assumes that the data‐generating process involves two linked stochastic processes: first a Bernoulli process (with or without covariates) describes the probability of observing a non‐zero; and subsequently a beta regression model is specified for the value of the proportion itself for all non‐zero observations. It can therefore be thought of as a special case of two‐component finite mixture modeling. This approach is also available for one‐augmented beta regression (Ospina & Ferrari, 2012) or zero‐and‐one augmented beta regression (Fang & Kong, 2015). See Joseph, Preston, and Johnson (2016) for an application of zero‐one‐augmented beta regression to vegetation cover, and Wright et al. (2017) for a further extension of zero‐augmented beta regression that leverages repeated observations to separately model true absences from apparent (observation‐error related) absences in vegetation survey data. A somewhat different method, involving modification of the likelihood function, has recently been proposed for zero‐augmented Dirichlet regression (Tsagris & Stewart, 2018).

The added value of a zero‐inflated/augmented models will be larger when both zeros and relatively high proportions are observed within replicates of the same treatment – evidence that two data‐generating processes are operating. Moreover, the separation of the data‐generating process into two independent components allows for additional inference about processes underlying presence/absence, as distinct from processes determining the observed proportions (see e.g. Keim, DeWitt, Fitzpatrick, & Jenni, 2017; Wright et al., 2017). Given these advantages, and the relative ease of working with zero‐augmented models within existing software packages (see Table 1), we suggest that zero/one‐augmented models should in general be used whenever there is a reasonable a priori expectation of zero or one values in the dataset, e.g. for species cover data. On the other hand, for certain data types (e.g. biomass partitioning in Case Study 2, below), zero or one values in the dataset will be absent in most cases.

Since the application of regular beta regression to data with zeros (and/or ones) requires transformation of the data, formal model selection criteria such as AIC or Bayesian Information Criterion (BIC) cannot be applied to compare the fit of a beta regression model fitted to a transformed response to zero‐and/or‐one inflated beta regression fit to an untransformed response. Therefore, model selection needs to be based on other criteria such as visual inspection of residuals and comparison of model predictions with observations. In case study 1, we compare the conclusions drawn from beta regression on transformed variables, and zero‐augmented beta regression (see also Appendix S3).

6.1 Case Study 1 Percent cover in quadrats

The first example involves experimental manipulation of the density of the sea urchin Centrosthepanus rodgersii to investigate its effect of grazing on the colonization of filamentous algae (Andrew & Underwood, 1993). Algae colonization was measured by percent cover in five 0.25 m² quadrats randomly located within larger patches subject to one of four levels of grazer removal treatment. Andrew and Underwood (1993) analysed this data (reanalysed by Quinn & Keough, 2002) using a nested ANOVA to account for the within‐patch replication. They concluded that treatments did not significantly affect the cover of filamentous algae.

In Appendix S3, we provide a detailed, step‐by‐step analysis of this dataset using different versions of beta regression, with accompanying code. We approach the analysis in two ways: first to illustrate the basic ideas we compare classical ANOVA to beta regression with and without a model for varying precision. To focus on the comparison of these basic model types, and the role of residuals in diagnostics, we use data pooled per patch and transform the data according to Equation (1) to initially avoid issues with nested observations, and the presence of zeroes, respectively. Secondly, we perform the analysis on the original data, retaining the hierarchical structure and comparing the results of models with and without zero‐augmentation. This approach is what we would recommend in a ‘real‐world’ analysis, since it incorporates a priori information about the presence of zeros and the structure of the experimental design.

Figure 2 and Table 2 compare the results of classical ANOVA (assuming normal errors), beta regression with a fixed ϕ, and beta regression with ϕ dependent on removal treatment. Model selection based on AIC clearly favours the variable ϕ model (Table 2). The improved fit is also evident from comparison of the residual plots (Figure 2, second column) – the first two models show a strong relationship between values of the standardized residuals and the fitted values. This is due to overestimation of the amount of variance in the control treatment plots, a fact also visible when comparing the posterior predictive densities of the first two models with the observed data (Figure 2, first column). Interestingly, there is not a large difference in the estimates for the mean of each group (Table 2), but the classical ANOVA model has much broader confidence intervals (leading to non‐significant pairwise comparisons, see Appendix S3) than the beta regression, and moreover predicts values outside the possible range of (0, 1). Pairwise comparisons of the groups based on the variable ϕ beta regression model indicated that the control treatment differed significantly from the other treatments, but the other treatments did not differ significantly from each other, both in terms of mean response and precision.

To illustrate the use of mixed‐effects and zero‐augmented beta regression, we analyse the original dataset in which replicate quadrats (observational units) are observed within each patch (experimental units). Given the non‐independence of quadrats within each patch a mixed‐effects model is fit with patch as the grouping variable and ϕ dependent on treatment. The predicted distributions for each treatment (Figure 3a) are broadly similar to those obtained from the variable ϕ model on the pooled data (Figure 2), however, the higher number of data points seems to have increased the precision of the estimates for the non‐control treatments. Incorporating a zero‐augmented component to the model, where the probability of observing a zero is modelled as a function of removal treatment, leads to only slightly adjusted posterior predictions for the mean of each group (Figure 3b,c). As expected, the main added value of the zero‐augmented model is a much more accurate prediction of zero observations in the dataset (Figure 3d). The inability of models without zero‐augmentation to reproduce this important feature of the dataset would limit their usefuless for making further predictions.

As for the analysis of the pooled data, the conclusions from both the mixed‐effects and zero‐augmented mixed effects models are that any form of any degree of sea urchin removal from this environment leads to an increase in algal cover. This finding contrasts with the conclusions in the original analyses of Andrew and Underwood (1993) and the reanalysis by Quinn and Keough (2002), both of which concluded that there was no significant difference in percentage cover of filamentous algae between treatments. We would argue that by choosing a more realistic model for the response variable, allowing the dispersion to vary between treatments, and explicitly modelling the occurence of zeroes, a beta regression model better captures the features of the dataset, and therefore provides a more reliable basis for inference.

6.2 Case study 2 Biomass partitioning in plants

The second dataset comes from a study testing whether differences in growth parameters between fast‐ and slow‐growing plant species at optimal nitrogen supply persisted at low nitrogen supply (Poorter & Sack, 2012; Poorter et al., 1995). Two species, Deschampsia flexuosa (slow growing) and Holcus lanatus (fast growing) were grown under low and high nitrate supply for a maximum of 49 days. Replicate individuals were harvested at regular intervals for determination of biomass in roots, stems and leaves. This case study is an illustration of proportions that arise in a situation where the size of the observational unit (total biomass) is not fixed, and where there are more than two continuous categories (biomass of stems, leaves and roots). Dirichlet regression (Appendix S4) was used, as the generalization of beta regression for situations where proportions are calculated for more than two categories.

The response variables were vectors of the proportions of total plant biomass in leaves (LMF), roots (RMF) and stems (SMF). These proportions were modelled as function of species identity and nitrate levels. In contrast to Poorter et al. (1995), we included time as a covariate to investigate the temporal dynamics of biomass partitioning. We refer to Poorter and Sack (2012) for other options regarding the analysis of biomass fractions. The most parsimonious model (based on AIC) explained the mean RMF and SMF as a function of time, a quadratic transformation of time, species, nitrate supply, total biomass, the interactions between species and nitrate supply, species and time, nitrate supply and time, and a three way interaction between species, nitrate supply and time. In addition, the precision was modelled as a function of species, nitrate supply, time and the interaction between species and time (Table 3).

How the different fixed effects determine LMF, SMF and RMF is difficult to infer from direct inspection of the parameter estimates of the best fitted model because the different fractions are interrelated. We therefore displayed the predicted values of this model in Figure 4. The predicted fractions show that in the high nitrate supply treatment, both species changed their allocation to shoots, roots and leafs in a similar fashion, while under low nitrate supply H. lanatus and D. flexuosa allocate root and shoot biomass differently over time. This explains the significant three way interaction between time, species and nitrate supply. Excluding the species term led to an increase in AIC of 412, emphasizing the importance of species–specific effects on biomass partitioning. No obvious pattern remained in Pearson residuals when plotted against the fitted values or the covariates.

Another way to gain insight into the species effect is to compute the ratio of organ biomass of the two species. This ratio can be calculated from the predicted biomass fractions of the Dirichlet regression model, and can be thought of as a measure of effect size expressing the relative partitioning of biomass invested in leaves, stems and roots in H. lanatus compared to D. flexuosa (Figure 5). For example, early on in development, it is predicted that H. lanatus invest up to 1.5 times more in roots than D. flexuosa, while at harvest the investment in roots is similar. Despite this pattern, the 95% prediction interval of the investment ratios, taking the variation of individual replicates into account, shows that the variation in investment ratio within species is substantial.

6.3 Case study 3 Percent cover and comparison of beta regression and transformations‐based approaches

In this case study, beta regression and transformation‐based approaches are compared to illustrate the mismatch between observations and predictions that can arise when using transformation‐based approaches or when choosing an inappropriate link function within beta regression. A synthetic dataset was used to compare the performance of various approaches against true ground cover.

We created a dataset of tree cover using a stochastic, two dimensional spatial model where tree density is modeled as a function of mean annual precipitation (mm/year) following findings of (Hirota, Holmgren, Nes, & Scheffer, 2011; Staver, Archibald, & Levin, 2011). Importantly, the underlying data‐generating process is not directly related to any of the model specifications we are comparing. Projected ground‐cover of a range of 20 forests was simulated that varied in the mean annual rainfall received (ranging from 125 to 2,500 mm/year). Trees were positioned randomly in the area by drawing coordinates from a uniform distribution within the grid, and the size of the (circular) individuals was simulated by sampling values of crown diameter from a lognormal distribution. Percent cover on 1 ha plots was ‘estimated’ by simulating 15 randomly positioned non‐overlapping quadrats of 10 × 10 m² (Figure 6a and Appendix S5 for details). We then used these simulated samples to estimate the relationship between mean percent cover and mean annual precipitation using one of five methods: log transformation or logit transformation followed by an ordinary least squares linear regression model, or a beta regression with either a cloglog, logit, or log link. To avoid fitting problems in the beta model with log link, we fitted a regression model on the proportion of non‐cover (i.e. 1 − cover) and set the intercept at 1. This constrains the model to return zero cover at values of zero annual precipitation, which is biologically plausible in this case.

The beta regression with log link best fitted the data (Figure 6). The root mean squared error (RMSE) of the model predictions versus true tree cover was on average 0.058, 0.073 and 0.099 for the log, logit and cloglog link respectively, and 0.069 and 0.663 for the logit and log transformation respectively (average of 50 simulations). In all cases, the RMSE decreased with increasing quadrat size while keeping the area sampled constant to 16% of the total area. However, the RMSE for the logit transformation decreased much faster compared to the beta regression models, but always stayed above the RMSE of the beta regression model with log link. The stronger reduction in RMSE with increasing quadrat size compared to the decrease in RMSE in beta regression models illustrates the point that larger variance in the observed proportions increases the mismatch between observations and the predictions of transformation‐based approaches (see Appendix S5). Furthermore, the choice of the link function substantially affected model fit (Figure 6). Thus, the simulations show that beta regression is better able to predict tree cover compared to transformation‐based approaches, provided that the link function is chosen carefully.